###############################################
#TABLES B1 - B3
###############################################

library(lmtest)
library(sandwich)
library(stargazer)

####load cleaned data####

data1 <- read.csv("./data/data1.csv")
data2 <- read.csv("./data/data2.csv")
data3 <- read.csv("./data/data3.csv")

####TABLE B1 - Identity Threat vs Control (data1)####

lm1a <- lm(prej_direct_1 ~ as.factor(trt) + education + income + newjob, data=data1)
lm1b <- lm(list_items ~ as.factor(trt) + education + income + newjob + as.factor(trt)*sens, data=data1)
lm1c <- lm(donation_amount ~ as.factor(trt) + education + income + newjob, data=data1)
lm1d <- lm(donation_dummy ~ as.factor(trt) + education + income + newjob, data=data1)
lm1e <- lm(natid_pca ~ as.factor(trt) + education + income + newjob, data=data1)
lm1f <- lm(empire_1 ~ as.factor(trt) + education + income + newjob, data=data1)
lm1g <- lm(comparison_1 ~ as.factor(trt) + education + income + newjob, data=data1)

lm1r <- coeftest(lm1a, vcov = vcovHC(lm1a, type = "HC0"))
lm2r <- coeftest(lm1b, vcov = vcovHC(lm1b, type = "HC0"))
lm3r <- coeftest(lm1c, vcov = vcovHC(lm1c, type = "HC0"))
lm4r <- coeftest(lm1d, vcov = vcovHC(lm1d, type = "HC0"))
lm5r <- coeftest(lm1e, vcov = vcovHC(lm1e, type = "HC0"))
lm6r <- coeftest(lm1f, vcov = vcovHC(lm1f, type = "HC0"))
lm7r <- coeftest(lm1g, vcov = vcovHC(lm1g, type = "HC0"))

robust_se1    <- lm1r[,2]
robust_se2    <- lm2r[,2]
robust_se3    <- lm3r[,2]
robust_se4    <- lm4r[,2]
robust_se5    <- lm5r[,2]
robust_se6    <- lm6r[,2]
robust_se7    <- lm7r[,2]

stargazer(lm1a, lm1b, lm1c, lm1d, lm1e, lm1f, lm1g, type = "latex",
          se = list(robust_se1, robust_se2, robust_se3, robust_se4, robust_se5, robust_se6, robust_se7),
          dep.var.labels = c("Pub. Prejudice", "List Items", "Donation", "Donation (Dummy)", "Nat. ID", "Empire", "Whataboutism"),
          covariate.labels = c("Treatment 1", "Education (Low)", "Income", "New Job", "Sensitive Item", "Sensitive Item * Treatment"), 
          model.numbers = F, column.sep.width = "0pt")

####TABLE B2 - Whataboutism vs Identity Threat (data2)####

lm2a <- lm(prej_direct_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm2b <- lm(list_items ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob + as.factor(trt)*sens, data=data2)
lm2c <- lm(donation_amount ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm2d <- lm(donation_dummy ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm2e <- lm(natid_pca ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm2f <- lm(empire_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)
lm2g <- lm(comparison_1 ~ as.factor(trt) + lr_1 + natid_covar1_1 + income + newjob, data=data2)

lm1r <- coeftest(lm2a, vcov = vcovHC(lm2a, type = "HC0"))
lm2r <- coeftest(lm2b, vcov = vcovHC(lm2b, type = "HC0"))
lm3r <- coeftest(lm2c, vcov = vcovHC(lm2c, type = "HC0"))
lm4r <- coeftest(lm2d, vcov = vcovHC(lm2d, type = "HC0"))
lm5r <- coeftest(lm2e, vcov = vcovHC(lm2e, type = "HC0"))
lm6r <- coeftest(lm2f, vcov = vcovHC(lm2f, type = "HC0"))
lm7r <- coeftest(lm2g, vcov = vcovHC(lm2g, type = "HC0"))

robust_se1    <- lm1r[,2]
robust_se2    <- lm2r[,2]
robust_se3    <- lm3r[,2]
robust_se4    <- lm4r[,2]
robust_se5    <- lm5r[,2]
robust_se6    <- lm6r[,2]
robust_se7    <- lm7r[,2]

stargazer(lm2a, lm2b, lm2c, lm2d, lm2e, lm2f, lm2g, type = "latex",
          se = list(robust_se1, robust_se2, robust_se3, robust_se4, robust_se5, robust_se6, robust_se7),
          dep.var.labels = c("Pub. Prejudice", "List Items", "Donation", "Donation (Dummy)", "Nat. ID", "Empire", "Whataboutism"),
          covariate.labels = c("Treatment 2", "L-R", "Nat. ID", "Income", "New Job", "Sensitive Item", "Sensitive Item * Treatment"), 
          model.numbers = F, column.sep.width = "0pt")

####TABLE B3 - Whataboutism vs Control (data3)####

lm3a <- lm(prej_direct_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3)
lm3b <- lm(list_items ~ as.factor(trt) + lr_1 + unemployed + as.factor(trt)*sens, data=data3)
lm3c <- lm(donation_amount ~ as.factor(trt) + lr_1 + unemployed, data=data3)
lm3d <- lm(donation_dummy ~ as.factor(trt) + lr_1 + unemployed, data=data3)
lm3e <- lm(natid_pca ~ as.factor(trt) + lr_1 + unemployed, data=data3)
lm3f <- lm(empire_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3)
lm3g <- lm(comparison_1 ~ as.factor(trt) + lr_1 + unemployed, data=data3)

lm1r <- coeftest(lm3a, vcov = vcovHC(lm3a, type = "HC0"))
lm2r <- coeftest(lm3b, vcov = vcovHC(lm3b, type = "HC0"))
lm3r <- coeftest(lm3c, vcov = vcovHC(lm3c, type = "HC0"))
lm4r <- coeftest(lm3d, vcov = vcovHC(lm3d, type = "HC0"))
lm5r <- coeftest(lm3e, vcov = vcovHC(lm3e, type = "HC0"))
lm6r <- coeftest(lm3f, vcov = vcovHC(lm3f, type = "HC0"))
lm7r <- coeftest(lm3g, vcov = vcovHC(lm3g, type = "HC0"))

robust_se1    <- lm1r[,2]
robust_se2    <- lm2r[,2]
robust_se3    <- lm3r[,2]
robust_se4    <- lm4r[,2]
robust_se5    <- lm5r[,2]
robust_se6    <- lm6r[,2]
robust_se7    <- lm7r[,2]

stargazer(lm3a, lm3b, lm3c, lm3d, lm3e, lm3f, lm3g, type = "latex",
          se = list(robust_se1, robust_se2, robust_se3, robust_se4, robust_se5, robust_se6, robust_se7),
          dep.var.labels = c("Pub. Prejudice", "List Items", "Donation", "Donation (Dummy)", "Nat. ID", "Empire", "Whataboutism"),
          covariate.labels = c("Treatment 2", "L-R", "Unemployed", "Sensitive Item", "Sensitive Item * Treatment"), 
          model.numbers = F, column.sep.width = "0pt")
